Monte Carlo simulations of bosonic reaction-diffusion systems 
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An efficient Monte Carlo simulation method for bosonic reaction-diffusion systems which are 
mainly used in the renormalization group (RG) study is proposed. Using this method, one- 
dimensional bosonic single species annihilation model is studied and, in turn, the results are com- 
pared with RG calculations. The numerical data are consistent with RG predictions. As a second 
application, a bosonic variant of the pair contact process with diffusion (PCPD) is simulated and 
shown to share the critical behavior with the PCPD. The invariance under the Galilean transfor- 
mation of this boson model is also checked and discussion about the invariance in conjunction with 
other models are in order. 

PACS numbers: 64.60.Ht, 05.10.Ln, 89.75.Da 
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I. INTRODUCTION 

The reaction-diffusion (RD) systems have become a 
paradigm for studying certain physical, chemical, and bi- 
ological systems [Q . In the study of the RD systems on 
a lattice via Monte Carlo (MC) simulations, particles in- 
volved in the dynamics usually have hard core exclusion 
property. In other words, MC simulations have been in- 
terested in the lattice systems where multiple occupancy 
at a lattice point is prohibited. These particles are often 
referred to as fermions, but this paper prefers the term 
"hard core particles." Meantime, the renormalization- 
group (RG) calculations that have been applied success- 
fully to several RD systems are in many cases based on 
the path integral formalism for classical particles with- 
out hard core exclusion, or, if we are allowed to abuse 
terminology, bosons On this account, the com- 

parison of the numerical studies to the RG calculations 
is sometimes nontrivial. 

There are two ways to fill a gap between numerical 
and analytical studies. One is to make a path integral 
formula for hard core particles which is suitable for the 
RG calculations. Actually this path has been sought 
and some formalisms are suggested 0, IE Q • The other 
is to find a numerical method to simulate boson sys- 
tems. In this context, numerical integration studies of 
the equivalent Langeyin equations to boson systems have 
been performed 0, LJ Hfl 0~D] • However, it is not always 
possible to find an equivalent Langevin equation [12j . 
By the same token, the applicability of this approach is 
somewhat restricted. Thus, another numerical method is 
called for. To our knowledge, no algorithm to simulate 
general bosonic RD systems directly has been suggested 
and to find such a algorithm is still a challenging topic. 

This paper suggests an algorithm to simulate the 
bosonic RD systems. Section |H] is devoted to a heuristic 
explanation of the algorithm to simulate general bosonic 
single species RD systems. In Sec. IIIII the numerical 



method applies to some bosonic RD systems. At first, 
the single species annihilation models with various con- 
ditions are simulated, along with the comparison to the 
RG predictions. Then, a bosonic version of the pair con- 
tact process with diffusion is discussed, focusing on the 
universality and Galilean invariance. Section [W] summa- 
rizes the works. 



II. ALGORITHM 

This section explains the algorithm suitable for MC 
simulations of bosonic RD systems. Although the discus- 
sion in this section is restricted to single species cases, the 
extension to multispecies problems is straightforward. 

The reaction dynamics of diffusing bosons is repre- 
sented as 



nA (n + m)A, 



(1) 



where n > 0, m > — n, m ^ 0, and X nm is the tran- 
sition rate. Each particle diffuses with rate D on a d 
dimensional hypercubic lattice. The periodic-boundary 
conditions are assumed, but other boundary conditions 
do not limit the validity of the algorithm. Configurations 
are specified by the occupation number p x (> 0) at each 
lattice point x. A configuration is denoted as {/?} which 
means {p^ja; £ L d }, where L d stands for the set of the 
lattice points and the cardinality of L d is L d . The master 
equation which describes stochastic processes modeled by 
Eq. (JTJ takes the form [II G3 
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where P = P({p},t) is the probability with which the 
configuration of the system is {p} at time t, (x,y) 
means the nearest-neighbor pair (x,y € L d ), f{p Xl n) — 
(Pxty/(Px — n)\ is the number of ordered n tuples at site 
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FIG. 1: An example of a configuration change of a one- 
dimensional RD system with L — 4 due to hopping. The black 
circle signifies a particle and numbers below the horizontal line 
indicate the lattice point. A particle at site 1 hops to site 2. 



x of the configuration {p}, and E xy and C Xjm are oper- 
ators affecting P({p},t) such that 



E x ^yP 



P = P({--- , Px -m,---};t). 



};*), 



(3) 



(4) 



The master equation implies that during infinitesimal 
time interval dt, the average number of transition events 
for the configuration {p} is 

E(dt,{p}) =dtJ2[ 2dD5 nA + J2 Km ) f(px, n) 

x.n \ m / 

= dty^i 2dDS n> i + ^n!A nm g(p x ,n), 

x,n \ m / 

where g(p X: n) = f(p x ,n)/n\ — ( p ^) is the number of 
(nonordered) n tuples at site x. Therefore, the first step 
for MC simulations is to select one of n tuples with an 
equal probability. For the convenience of description and 
better understanding, we introduce a model dependent 
function h(p x ,n) = e n g(Px,n), where e„ takes 1 (0) if 
D§ n> i +J2 m ^ nm ^ s nonzero (zero). The meaning of e n is 
straightforward; we do not have to consider the reaction 
dynamics with transition rate zero (see below). 

The simplest way to implement the selection is as fol- 
lows: First a site x is selected with probability N x /M, 
where N x — ^2 n h(p x , n) which will be called the number 
of accessible states at site x and M = J2 X N x ■ Then, 
n is chosen with probability h(p x , n)/N x which is zero if 
e„ = 0. For this procedure, the array of the number of 
particles at all sites, say p[ ] (p[x] = p x ), is necessary. 

However, it is not efficient because there are too many 
floating number calculations. For a faster performance, 
we introduce two more arrays, say list[ ] and act[ ][ ]. 
The array list[ ] refers the location of any n tuple. Each 
element of list[ ] takes the form (x,£), where s is a 
site index and £ lies between 1 and N x (1 < £ < N x ). 
From £ and p[x], which n tuple is referred by the array 
list[ ] is determined. If £ < h(p x ,0), n = is implied. 
Else if £ < h(p x ,0) + h(p x , 1), n = 1 is meant. Else if 
I < h(p X7 0) + h(p x , 1) + h(p x , 2), I indicates one of pairs 
at site x, and so on. In case the total number of accessible 
states in the system is M, the size of list[ ] is M and all 
elements of list[ ] should satisfy that list[p] 7^ list[g] 
if p ^ q (1 < p,q < M). Hence, the random selec- 
tion of an integer between 1 and M is equivalent to the 



choice of one n tuple among M accessible states with 
an equal probability. The array act [ ] ] is the inverse 
of the list[ ], that is, list[s] = (x,£) corresponds to 
act[cc][f] = s. 

After selecting x and n, the transition nA — > (n + m)A 
occurs with the probability of n\X nm At for all possible m, 
where At is independent from configurations. Provided 
n = 1 is selected, in addition to reaction processes, a par- 
ticle at x hops to one of the nearest neighbors with prob- 
ability DAt. To make the transition probability have a 
meaning, At should satisfy 



2dDS nA 



At < 1, 



(5) 



for all n. Time is increased by At/M . On average, this 
algorithm generates E(At, {p}) transition events during 
time interval At. After the system's evolving, three ar- 
rays, p, list, and act, are updated in a suitable way 
(see below). 

Through an example, how the system evolves in silico 
is to be clarified. Consider a RD system with e n = 
for n > 3 and n — 0. In this case, N x — Y] h(p x ,n) — 
g(p x , 1) -I- g(p x ,2) = p x (p x + l)/2 will be used. Assume 
that we are given a configuration p[l] = 2, p[2] = 0, 
p[3] = 1, and p[4] = 3 (JVi = 3, N 2 = 0, = 1, - 6, 
hence M = 10); see Fig. Q] Complete lists of two arrays 
list[ ] and act[ ][ ] for this configuration are illustrated 
on the left-hand side of Table [I] The algorithm starts 
from selecting one number between 1 and M, randomly. 
Let us assume that 2 is selected, which makes list [2] 
to be checked. Since list[2] = (1,2) and 2 < p[l], a 
particle dynamics at site 1 will be attempted. Again as- 
sume that a hopping to the site 2 whose probability is 
DAt occurs, which results in a change of the configura- 
tion as shown in Fig. ^ Accordingly, three arrays should 
be updated. Figure [21 shows how the evolution is coded 
(based on the language c). In this code, rho[a;] is the 
number of particles at site x (— p x ), N[x] is the number 



TABLE I: An example of making two arrays referring each 
other from the configuration shown in Fig. Q Two columns 
on the left (right) hand side correspond to the configuration 
before (after) the hopping event. 



Before 



After 



list[l] = 


(1,1) 


act[l][l] 


= 1 


list[l] 


= (1,1) 


act[l][L> 


= 1 


list [2] = 


(1,2) 


act[lp; 


=2 


list [2] 


= (4,6) 


act[2][L> 


=9 


list[3] = 


(1,3) 


act[l][3] 


=3 


list [3] 


= (4,5) 


act[3][l] = 


=4 


list[4] = 


(3,1) 


act[3][l 


=4 


list [4] 


= (3,1) 


act[4][l]= 


=5 


list [5] = 


(4,1) 


act[4][l] 


=5 


list [5] 


= (4,1) 


act [4] [2] = 


=6 


list[6] = 


(4,2) 


act [4] [2 


=6 


list [6] 


= (4,2) 


act [4] [3] = 


=7 


list [7] = 


(4,3) 


act [4] [3 


=7 


list [7] 


= (4,3) 


act [4] [4] = 


=8 


list [8] = 


(4,4) 


act [4] [4 


=8 


list [8] 


= (4,4) 


act [4] [5] = 


=3 


list[9] = 


(4,5) 


act [4] [5 


=9 


list [9] 


= (2,1) 


act [4] [6] = 


=2 


list[10] 


= (4,6) 


act [4] [6] 


=10 
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f or (k=N [rho [1] -1] +1 ;k<=N [rho [1] ] ; k++) { 
s = act [1] [k] ; 

x = list [M] [0] ; L = list[M] [1] ; 

act [x] [L] =s ; 

act[l] [k]=0; 

list [s] [0]=list [M] [0] ; 

list [s] [l]=list[M] [1] ; 

M=M-1; 

> 

rho[l]=rho[l]-l; 

for (k=N [rho [2] ] +1 ; k<=N [rho [2] +1] ; k++) { 
M=M+1 ; 

act [2] [k]=M; 

list[M] [0]=2;list[M] [l]=k; 

> 

rho[2]=rho[2]+l; 

FIG. 2: A program which updates three arrays p\ ], list[ ], 
and act [ ] [ ] after the hopping event shown in Fig. Q 



of accessible states at site x (= N x ), and each element 
of list [ ] is treated as an array. The first (second) for 
loop signifies the decreasing (increasing) of the number 
of accessible states at site 1 (2), which can be used for 
any particle number decreasing (increasing) events. The 
code generates the lists on the right-hand side of Table 
[I] Time is increased by At/ 10. Then again choose one 
number between 1 to 9, randomly, and so on. 

Equipped with the numerical methods, Sec. Illll studics 
some bosonic RD systems which show scaling behavior. 



III. APPLICATIONS 

A. Single species annihilation model 

The algorithm explained in the previous section is 
applied to a one-dimensional single species annihilation 
model which corresponds to X nm — unless n = 2 and 
m = —2. For saving the writing effort, let us rename 
A2,-2 A. The renormalization-group calculation pre- 
dicts that the annihilation fixed point corresponds to 
A = oo ■ Infinite pair annihilation rate means that two 
particles occupying the same site by any chance will be 
removed instantaneously. Accordingly, at most one parti- 
cle can reside at each site. Hence, the boson model with 
infinite annihilation rate is equivalent to the diffusion- 
limited annihilation model (DLAn) of hard core particles 
which can be solved exactly 0]. It * s known that the 
particle density of the DLAn starting from the random 
initial condition decays as 



Pit) 



lim 



1 L 

7 



l 



x=l 



(1 + 0(1/*)). (6) 



This behavior does not depend on the initial density. 
Since renormalized coupling constant flows to the annihi- 
lation fixed point, the asymptotic behavior of the density 




FIG. 3: (Color online) A log-log plot of p(t)V8TvDt as a 
function of t for various D with A = | and po = 1. All curves 
approach to 1 as t increases. Inset: same, but density is not 
multiplied by \fSixDt. 



for finite A is expected to be the same as Eq. Besides, 
it is expected that the smaller the value of A, the later the 
system enters the scaling regime. Actually, these predic- 
tions are tested for the annihilation model of hard core 
particles 0]. However, to our knowledge, there is no 
satisfactory numerical test for the RG predictions using 
a boson model [l6| . 

The Poisson distribution is used as an initial condition, 
which can be implemented if we randomly distributed 
PqL particles on the lattice. For this distribution, the 
probability that q particles reside at site x is 



Px(q) = 




poL—q 



4 e ~ p °' ( 7 ) 



where L is assumed to be sufficiently large and q <C poL. 
Using the algorithm explained in the previous section and 
varying D, A, and po, we simulated the one-dimensional 
annihilation model. The system size is 2 16 and the num- 
ber of independent samples is 200 for each data set. 
Figure[3]shows the decaying behavior of the density for 
- , ~, and jt- with p = 1 and A = i. Each curves 
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FIG. 4: (Color online) A log-log plot of p{t)y/Aiit vs t for 
various po with D — i and A = i. In the asymptotic regime, 
all data sets show the same behavior. 
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FIG. 5: (Color online) A log-log plot of p(t)\/A-Kt vs t for 
various A with D = | and po = 1. Although the system 
with smaller A enters the scaling regime slowly, all curves 
eventually meet for large t. 



approaches to l/y/SnDt as the RG calculation predicted. 
We also check the initial condition dependence, by sim- 
ulating systems with various initial density 2, 1, i, and 
2 with D = h and A = h . Figure 01 shows the initial 
condition independence of the asymptotic behavior. Fi- 
nally, we also confirm that the asymptotic behavior is 
not affected by A, see Fig. GD As expected, the system 
with smaller A enters the scaling regime later. The MC 
simulation for bosonic annihilation models confirms the 
predictions of the RG study 



B. Pair contact process with diffusion 

The pair contact process with diffusion (PCPD) is a 
RD system of diffusing hard core particles with two com- 
peting dynamics of 2A — > 3A (fission) and 2A — > (an- 
nihilation) , which shows a continuous transition ■ At 
first sight, the bosonic variant of the PCPD might be 
regarded as the boson model with X nm = except A21 
and A2.-2- However, this variant does not show a contin- 
uous transition and there is no steady state in its active 
(fission dominating) phase To have a well-defined 

steady state in all phases, a mechanism to keep the den- 
sity from blowing up is required. Introducing a triple 
reaction such as 3A — > 2A, one can get a model with 
well-defined steady states. Although the boson model 
with \ nm = except A21, \2.-2, and A3.-1 has been 
expected to show a continuous transition [T3 , MC simu- 
lation results for this type of boson model which will be 
called "BPCPD" has yet been reported in the literature, 
although a parallel update bosonic model with so-called 
soft-constraint was studied [Tsj . 

Using parameter values D = ~, X^^-i = g, \2,-2 = 
p/2, and A21 = (l—p)/2 where p is the tuning parameter, 
the critical behavior of the BPCPD is studied. As an 
initial condition, we set p x — 2 for all x (1 < x < L). 
Figure E] shows the decaying behavior at criticality of two 
order parameters, the particle and pair densities which 
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FIG. 6: (Color online) Time dependence of the particle 
and pair densities multiplied by t' 3 '"!! with f3/v\\ = 0.205 in 
semilogarithmic plot at criticality for the BPCPD. Inset: ef- 
fective exponents of the order parameters at p — 0.148 79. 



are defined as 



Pi 



Pi 



X 

(*) = 7 X>« (*.-!)>« 



(8) 



where (•••)* means the average over ensembles at time 
t. The system size in use is 2 15 and all samples (around 
10 3 samples are independently simulated) up to obser- 
vation time (~ 2.5 x 10 6 MC steps) have at least one 
site with two or more particles. The critical point is 
found to be p c = 0.148 79(1) with the critical exponent 
ft/vn = 0.205(5) which is estimated from the effective 
exponent 



-S(t) 



ln[pi, 2 (t)] -ln[pi, 2 (t/m)] 
lnm 



(9) 



with m = 10. At criticality, S(t) approaches to ft/vn as 
t goes to infinity. The simulation results are consistent 
with the previous works within error bars [p^ . l2fl |. Hence, 
we conclude that the BPCPD has the same critical scal- 
ing with the PCPD. 

Following the path integral formalism for bosonic RD 
systems pj, the action of the BPCPD, S = J dt dx£, 
after taking the (naive) space-time continuum limit has 
the form 



C = (/)(dt-DV 2 



-9i # +A 3 <; 



(10) 



which is the same as one studied in Ref. |2l| which is 
derived from path integral formalism for the exclusive 
particle systems introduced in Ref. 0. It is argued, 
however, via RG calculations [2l| and numerical studies 
[20I l2l| that Eq. (|10|) is inappropriate for studying the 
critical behavior of the PCPD using the RG techniques. 
Nonetheless, we will show that the Galilean invariance 
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FIG. 7: (Color online) Space-time configuration for the unbi- 
ased and biased BPCPD models at criticality. Blue (red) dots 
represent the sites with only one particle (at least two parti- 
cles) and white dots stand for the empty sites, (a) and (b) 
are configurations of the BPCPD with Dr — i and Dr = 1, 
respectively, (c) is the same as (b) except that the space 
coordinate is Galilean transformed with v — Dr — Dl — 1. 



FIG. 8: (Color online) Space-time configuration for (a) the 
PCPD and (b) the DPCPD at criticality studied in Ref. 
prU . Blue (red) dots represent the isolated particles (par- 
ticles which are members of pairs) and white dots stand for 
the empty sites, (c) is the same as (b) except that the space 
coordinate is Galilean transformed with v — 1. 



(GI) of the BPCPD, which is anticipated from Eq. JTDJ, 
is still correct in the strong sense (see below). 

For some RD systems, biased diffusion only changes 
nonuniversal constants such as the critical point and does 
not affect the critical behavior. Examples are the driven 
branching annihilating random walks (DBAW) studied 
in Ref. [20|. Such systems will be called to have the GI 
in the weak sense (Glweak). Why the critical point is 
dependent on the bias strength is understandable within 
the framework of Ref. pj. Using the path integral for- 
malism for hard core particles introduced in Ref. |?J , the 
terms appearing in the action due to the bias with the 
strength v take the form 

£bias = V (<&e0||0a> - 4i<f>xd\\<(>as) , (11) 

where du is the lattice gradient defined as dncf> x = 
{<p<c+e« ~ ( / , x-e||)/2 with e || the unit vector along the 
bias direction. The derivation of Eq. (|llfl is shown in 
the Appendix. The Galilean transformation gauges away 
the first term in Eq. 111(1. but cannot remove the second 
term. Since the second term in Eq. I jlljl is irrelevant in 
the RG sense for the DBAW, this does not affect the 
universal behavior, but the very existence of this irrele- 
vant term can change the critical point. Therefore, the 
DBAW is of the Glweak. Meanwhile, the PCPD is not 
of the GI even in the weak sense [2(| • Since it is shown 
that the field theory with the action (|10|) is not viable 
|2l|. we cannot extract any information from Eq. i|ll|) 
concerning the driven pair contact process with diffusion 
(DPCPD). To understand the DPCPD and the PCPD 
from the field theoretical point of view, more elaborated 
studies are required. 

The bias diffusion of bosons does not generate the sec- 
ond term in Eq. (|llf> . In this context, the Galilean trans- 
formation totally gets rid of the effect of bias for bosons. 
Hence, two systems with or without bias have the same 
probability distribution, let alone the critical behavior. 
These systems will be called to have the GI in the strong 
sense (GIstrong). Consider a one-dimensional bosonic 



RD system with reaction dynamics in Eq. Q in which 
each particle hops to the right (left) with rate Dr (Dl). 
The GIstrong for this model means that whatever value 
Dr takes with the constraint Dr + Dl = 2D (constant), 
the system shares the probability distribution with the 
unbiased model (Dr — Dl = D). It is checked numer- 
ically for various Dr with Dr + Dl = 1, whether the 
BPCPD has the GIstrong or not. We observed that the 
particle and pair densities have the same behavior at the 
same p within statistical error (not shown here) . In Fig. 
space-time configurations of the BPCPD models with 
the unbiased diffusion (Dr = D = |), fully biased dif- 
fusion (Dr = 2D = 1), and the Galilean transformation 
for the full bias case are shown. After Galilean trans- 
formation, no noticeable difference between biased and 
unbiased cases is observed. For comparison, we present 
in Fig. [SI the space-time configuration of the PCPD and 
the DPCPD studied in Ref. 20]. As the Galilean trans- 
formed space-time configuration shows, the bias cannot 
be removed in the DPCPD. The Galilean transformation 
generates the biased motion of the paired particles which 
shows the existence of the relative bias between isolated 
particles and paired ones. Although the validity of Eq. 
((10JI as an appropriate action for the RG study regard- 
ing the PCPD is rather problematic, any single species 
bosonic RD systems with on-site reactions are conjec- 
tured to have the GIstrong. 

The discussion about the GIstrong should be restricted 
to boson models with random sequential update dynam- 
ics. If the dynamics occurs in a parallel way as in Ref. 
[T^|. the GI argument from the invariance of the local 
action like Eq. (|10f> under the Galilean transformation is 
not directly applicable. Even worse, the one dimensional 
system with pr = 1 (for the definition of pr, see the next 
paragraph) is reduced to a single-site problem which is 
not expected to show phase transition. Notwithstanding, 
except this pathological case, the soft-constraint PCPD 
(SCPCPD) studied in Ref. is expected to have the 
Glweak [23. 

To understand what is happening in the SCPCPD, let 
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us explain the dynamics of the model. During unit time, 
changes of a configuration occur in two steps. At first, 
every particle hops to the right (left) with probability pr 
(j>l) and stays still with probability^ (pr+Pl+Ps = 1)- 
In Ref. 0], pr = pl = \ and ps — are used. After 
the hopping events, reactions occur at all sites. Rather 
interestingly, the model with pl = is statistically equiv- 
alent to the system with ps = provided pr is the same. 
When ps — 0, particles at the even sites do not inter- 
act with those in the odd sites. For example, see Fig. 
and regard the left figure of it as a configuration for the 
SCPCPD with ps — under the condition of the periodic 
boundary. At the end of the hopping event, particles at 
sites 1 and 3 (2 and 4) move on to sites 2 and 4 (1 and 
3). Thus, a system with size 2L (let us call it system A) 
can be considered two independent systems with size L 
(call it system B), if we interpret the hopping events to 
the left in the system A as a staying event in the system 
B. Since the system with pl = has a bias effect in 
diffusion except the pathological case of pr — 1, the GI 
for the SCPCPD is in a sense predictable. 

As a final remark, we would like to mention how the 
DPCPD behavior can be observed in the BPCPD model. 
As explained before, the bias applied to all particles has 
no effect. As was done for the SCPCPD in Ref. [U, 
if different bias is applied to a particle at singly occu- 
pied sites and a particle at multiply occupied sites, the 
DPCPD behavior such as mean-field-like exponents, log- 
arithmic corrections, etc., was observed (not shown here). 
This unusual bias cannot be included in the action like 
Eq. I|1U|) in a simple way, so this DPCPD behavior is not 
contradictory to the GIstrong of the BPCPD. 



IV. 



SUMMARY 



To summarize, an efficient algorithm is proposed to 
simulate the general bosonic reaction-diffusion systems 
and applies to the single species annihilation model and 
the bosonic variant of the pair contact process with dif- 
fusion. For the single species annihilation model, renor- 
malization group predictions are confirmed numerically. 
The BPCPD model is found to belong to the PCPD uni- 
versality class and maintains the Galilean invariance in 
the strong sense. Due to the lack of the analytical predic- 
tions for the PCPD, only the comparison of our results 
to published simulation results are possible. 
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APPENDIX: DERIVATION OF Eq. itTTl) 

From the path integral formalism for RD systems of 
hard core particles introduced in Ref. 0, Eq. will 
be derived in this appendix. Since the master equation 
is linear and the formalism in does not mix different 
dynamics, it is enough to consider the diffusion of hard 
core particles. For more detailed accounts, see Ref. 0. 

In general, the master equation becomes the imagi- 
nary time Schrodinger equation with (in general non- 
Hermitian) Hamiltonian H such that 



^\P;t) = -H\P;t), 



(A.I) 



where \P;t) = £ M P({p},t)\{p}) and \{p}) = U x \p x ) 
with p x taking cither 1 (occupied) or (vacant). To write 
down the Hamiltonian, introduced are the creation and 
annihilation operators for hard core particles in single 
species models which satisfy the following commutation 
relations: 

{&l,a x } = l, {a x ,a x } = {a x ,a x } = 0, 

[a a ,a a >] = [a x ,a x ,] = 0. 

Actually, these operators are nothing but the Pauli ma- 
trices. Using creation and annihilation operators, terms 
appearing in the Hamiltonian due to diffusion of hard 
core particles in the single species RD systems can be 
written as Hd = ^2 X H x with 



H -r. — 



d 

n 

i=l 



D 



2 



D-& 



v 



(n x v x+ei - a x a J a 

(fl x V x — ei o, x a x _ 



(A.3) 



the number operator, v x = 1 — n x , d 
is the unit vector along i direction, and hopping is biased 
along the || direction. 

The differential equation of the generating function F 
which is defined as 

F(mt) = J2 [ n*£- J p {{p\;t) = m\P;t), 



where 



takes the form 



-F = -{{(p}\H\P-t). 



(A.4) 



(A.5) 



(A.6) 



The generating function (|A.4|) corresponds to Eq. (15) 
of Ref. with the prescription (18a) in Ref. 0- Since 

({^}|«x = 6b (1 - <A^x)(MI, (A.7) 

(&}\a x = $ x {{Cp}\, (A.8) 

= <Px<p x {{<p}\, (A.9) 

{{(p}\v x = {l-Cp x (p x ){{ip}\, (A.10) 
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where (p-e = d/dtp x , one can find the partial differential 
equations for the generating function such that 



dt 



F = -£(M, {<f})F, 



(A.ll) 



with normal ordered evolution operator C which reads 



where V£ is the lattice Laplacian defined as V|/(x) — 
^2i = i(f (x + e{) + f (x — ei) — 2f(x)), and 3|| is the lattice 
gradient along the || direction defined as d\\f(x) — (f(x + 
en) — f(x — eii))/2. Since Eq. IjA.lip is a linear equation, 
we can write down the path integral solution with the 
action Q 



£(M, {<£}) = \ v [v^WP* ~ f%>'Px d \\ < P*] 



D 



+ (terms due to reactions), 



S= / &[$%</> + £({$},{</>})], 



(A.13) 



(A. 12) which completes the derivation of Eq. Qll[). 
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